mutate(model = data %>% map(~ lmer(income_equality ~ 1 + (1 | year) + (1 | iso3n) ,
data = .,
na.action = na.omit,
control=lmerControl(optCtrl=list(maxfun=6000)))),
tidied = model %>% map(~ tidy(., conf.int = TRUE)),
glance = model %>% map(~ glance(.)))
v1
#Inspect the data for imputation dataset 1
v1 %>%
filter(m == "imp1") %>%
unnest(tidied)
# Create a wide data frame of just the coefficients and standard errors
v1_par <- v1 %>%
unnest(tidied) %>%
select(m, term, estimate, std.error) %>%
gather(key, value, estimate, std.error) %>%
group_by(term, m) %>%
mutate(obs = row_number()) %>%
spread(term, value)
v1_par
v1_par <- v1_par %>%
drop_na(`(Intercept)`) %>%
select(-obs)
# Extract the coefficients
v1_coefs <- v1_par %>%
filter(key == "estimate") %>%
ungroup() %>%
select(-m, -key)
v1_coefs
# Extract the standard errors
v1_ses <- v1_par %>%
filter(key == "std.error") %>%
ungroup() %>%
select(-m, -key)
#Mi meld function Amelia
v1_melded <- mi.meld(v1_coefs, v1_ses)
v1_melded
#Regression coefficient table
v1_model_dg <- v1 %>%
unnest(glance) %>%
filter(m == "imp1") %>%
pull(df.residual)
v1_table <- as.data.frame(cbind(t(v1_melded$q.mi),
t(v1_melded$se.mi))) %>%
magrittr::set_colnames(c("estimate", "std.error")) %>%
mutate(term = rownames(.)) %>%
select(term, everything()) %>%
mutate(statistic = estimate / std.error,
conf.low = estimate + std.error * qt(0.025, v1_model_dg),
conf.high = estimate + std.error * qt(0.975, v1_model_dg),
p.value = 2 * pt(abs(statistic), v1_model_dg, lower.tail = FALSE))
v1_table
## Model 2
v2 <- all_imputations %>%
mutate(model = data %>% map(~ lmer(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
(1 | year) + (1 | iso3n) ,
data = .,
na.action = na.omit,
control=lmerControl(optCtrl=list(maxfun=6000)))),
tidied = model %>% map(~ tidy(., conf.int = TRUE)),
glance = model %>% map(~ glance(.)))
v2
#Inspect the data for imputation dataset 1
v2 %>%
filter(m == "imp1") %>%
unnest(tidied)
# Create a wide data frame of just the coefficients and standard errors
v2_par <- v2 %>%
unnest(tidied) %>%
select(m, term, estimate, std.error) %>%
gather(key, value, estimate, std.error) %>%
group_by(term, m) %>%
mutate(obs = row_number()) %>%
spread(term, value)
v2_par
v2_par <- v2_par %>%
drop_na(`age`) %>%
select(-obs)
# Extract the coefficients
v2_coefs <- v2_par %>%
filter(key == "estimate") %>%
ungroup() %>%
select(-m, -key)
v2_coefs
# Extract the standard errors
v2_ses <- v2_par %>%
filter(key == "std.error") %>%
ungroup() %>%
select(-m, -key)
#Mi meld function Amelia
v2_melded <- mi.meld(v2_coefs, v2_ses)
v2_melded
#Regression coefficient table
v2_model_dg <- v2 %>%
unnest(glance) %>%
filter(m == "imp1") %>%
pull(df.residual)
v2_table <- as.data.frame(cbind(t(v2_melded$q.mi),
t(v2_melded$se.mi))) %>%
magrittr::set_colnames(c("estimate", "std.error")) %>%
mutate(term = rownames(.)) %>%
select(term, everything()) %>%
mutate(statistic = estimate / std.error,
conf.low = estimate + std.error * qt(0.025, v2_model_dg),
conf.high = estimate + std.error * qt(0.975, v2_model_dg),
p.value = 2 * pt(abs(statistic), v2_model_dg, lower.tail = FALSE))
v2_table
## Model 3 ##
v3 <- all_imputations %>%
mutate(model = data %>% map(~ lmer(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
autocracy_5 + v2pepwrsoc_5 + autocracy_5*v2pepwrsoc_5 + v2x_polyarchy +
e_migdppcln + v2pepwrsoc +
(1 | year) + (1 + cohortmatch5_15 | iso3n) ,
data = .,
na.action = na.omit,
control=lmerControl(optCtrl=list(maxfun=6000)))),
tidied = model %>% map(~ tidy(., conf.int = TRUE)),
glance = model %>% map(~ glance(.)))
v3
#Inspect the data for imputation dataset 1
v3 %>%
filter(m == "imp1") %>%
unnest(tidied)
# Create a wide data frame of just the coefficients and standard errors
v3_par <- v3 %>%
unnest(tidied) %>%
select(m, term, estimate, std.error) %>%
gather(key, value, estimate, std.error) %>%
group_by(term, m) %>%
mutate(obs = row_number()) %>%
spread(term, value)
v3_par
v3_par <- v3_par %>%
drop_na(`ageCWC`) %>%
select(-obs)
# Extract the coefficients
v3_coefs <- v3_par %>%
filter(key == "estimate") %>%
ungroup() %>%
select(-m, -key)
v3_coefs
# Extract the standard errors
v3_ses <- v3_par %>%
filter(key == "std.error") %>%
ungroup() %>%
select(-m, -key)
#Mi meld function Amelia
v3_melded <- mi.meld(v3_coefs, v3_ses)
v3_melded
#Regression coefficient table
v3_model_dg <- v3 %>%
unnest(glance) %>%
filter(m == "imp1") %>%
pull(df.residual)
v3_table <- as.data.frame(cbind(t(v3_melded$q.mi),
t(v3_melded$se.mi))) %>%
magrittr::set_colnames(c("estimate", "std.error")) %>%
mutate(term = rownames(.)) %>%
select(term, everything()) %>%
mutate(statistic = estimate / std.error,
conf.low = estimate + std.error * qt(0.025, v3_model_dg),
conf.high = estimate + std.error * qt(0.975, v3_model_dg),
p.value = 2 * pt(abs(statistic), v3_model_dg, lower.tail = FALSE))
v3_table
v3_data <- v3_table %>%
drop_na()
v3_data
## Model 4 ##
v4 <- all_imputations %>%
mutate(model = data %>% map(~ lmer(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
autocracy_5 + v2dlencmps_5 + autocracy_5*v2dlencmps_5 + v2x_polyarchy +
e_migdppcln + v2dlencmps +
(1 | year) + (1 + cohortmatch5_15 | iso3n) ,
data = .,
na.action = na.omit,
control=lmerControl(optCtrl=list(maxfun=6000)))),
tidied = model %>% map(~ tidy(., conf.int = TRUE)),
glance = model %>% map(~ glance(.)))
v4
#Inspect the data for imputation dataset 1
v4 %>%
filter(m == "imp1") %>%
unnest(tidied)
# Create a wide data frame of just the coefficients and standard errors
v4_par <- v4 %>%
unnest(tidied) %>%
select(m, term, estimate, std.error) %>%
gather(key, value, estimate, std.error) %>%
group_by(term, m) %>%
mutate(obs = row_number()) %>%
spread(term, value)
v4_par
v4_par <- v4_par %>%
drop_na(`ageCWC`) %>%
select(-obs)
# Extract the coefficients
v4_coefs <- v4_par %>%
filter(key == "estimate") %>%
ungroup() %>%
select(-m, -key)
v4_coefs
# Extract the standard errors
v4_ses <- v4_par %>%
filter(key == "std.error") %>%
ungroup() %>%
select(-m, -key)
v4_ses
#Mi meld function Amelia
v4_melded <- mi.meld(v4_coefs, v4_ses)
v4_melded
#Regression coefficient table
v4_model_dg <- v4 %>%
unnest(glance) %>%
filter(m == "imp1") %>%
pull(df.residual)
v4_table <- as.data.frame(cbind(t(v4_melded$q.mi),
t(v4_melded$se.mi))) %>%
magrittr::set_colnames(c("estimate", "std.error")) %>%
mutate(term = rownames(.)) %>%
select(term, everything()) %>%
mutate(statistic = estimate / std.error,
conf.low = estimate + std.error * qt(0.025, v4_model_dg),
conf.high = estimate + std.error * qt(0.975, v4_model_dg),
p.value = 2 * pt(abs(statistic), v4_model_dg, lower.tail = FALSE))
v4_table
v4_data <- v4_table %>%
drop_na()
v4_data
## Model 5 ##
v5 <- all_imputations %>%
mutate(model = data %>% map(~ lmer(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
autocracy_5 + v2dlencmps_5 + v2pepwrsoc_5 +
autocracy_5*v2dlencmps_5*v2pepwrsoc_5 +
v2x_polyarchy +
e_migdppcln + v2dlencmps + v2pepwrsoc +
(1 | year) + (1 + cohortmatch5_15 | iso3n) ,
data = .,
na.action = na.omit,
control=lmerControl(optCtrl=list(maxfun=6000)))),
tidied = model %>% map(~ tidy(., conf.int = TRUE)),
glance = model %>% map(~ glance(.)))
v5
#Inspect the data for imputation dataset 1
v5 %>%
filter(m == "imp1") %>%
unnest(tidied)
# Create a wide data frame of just the coefficients and standard errors
v5_par <- v5 %>%
unnest(tidied) %>%
select(m, term, estimate, std.error) %>%
gather(key, value, estimate, std.error) %>%
group_by(term, m) %>%
mutate(obs = row_number()) %>%
spread(term, value)
v5_par
v5_par <- v5_par %>%
drop_na(`ageCWC`) %>%
select(-obs)
# Extract the coefficients
v5_coefs <- v5_par %>%
filter(key == "estimate") %>%
ungroup() %>%
select(-m, -key)
v5_coefs
# Extract the standard errors
v5_ses <- v5_par %>%
filter(key == "std.error") %>%
ungroup() %>%
select(-m, -key)
v5_ses
#Mi meld function Amelia
v5_melded <- mi.meld(v5_coefs, v5_ses)
v5_melded
#Regression coefficient table
v5_model_dg <- v5 %>%
unnest(glance) %>%
filter(m == "imp1") %>%
pull(df.residual)
v5_table <- as.data.frame(cbind(t(v5_melded$q.mi),
t(v5_melded$se.mi))) %>%
magrittr::set_colnames(c("estimate", "std.error")) %>%
mutate(term = rownames(.)) %>%
select(term, everything()) %>%
mutate(statistic = estimate / std.error,
conf.low = estimate + std.error * qt(0.025, v5_model_dg),
conf.high = estimate + std.error * qt(0.975, v5_model_dg),
p.value = 2 * pt(abs(statistic), v5_model_dg, lower.tail = FALSE))
v5_table
####Extract Texreg Table ####
extract_broom <- function(tidy_model, glance_model) {
# get estimates/standard errors from tidy
coef <- tidy_model$estimate
coef.names <- as.character(tidy_model$term)
se <- tidy_model$std.error
pvalues <- tidy_model$p.value
tr_object <- texreg::createTexreg(coef.names = coef.names,
coef = coef,
se = se,
pvalues = pvalues)
return(tr_object)
}
v1_tex <- extract_broom(v1_table)
v2_tex <- extract_broom(v2_table)
v3_tex <- extract_broom(v3_table)
v4_tex <- extract_broom(v4_table)
v5_tex <- extract_broom(v5_table)
library(texreg)
texreg(list(v1_tex, v2_tex, v3_tex, v4_tex, v5_tex),
head.tag = TRUE, body.tag = TRUE,
digits = 2,
booktabs = TRUE,
use.packages = FALSE,
caption = "Predicting Redistributive Preferences (Multiple Imputation Dataset)",
fontsize = "scriptsize",
label = "tab:tab-01")
View(all_imputations)
all_imputations$imp1
View(a.wvsdata)
data_imp1 <- a.wvsdata$imputations[[1]]
Sys.time()
ptm <- proc.time() # Start the clock
m3.ols <- brm(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
autocracy_5 + v2pepwrsoc_5 + autocracy_5*v2pepwrsoc_5 + v2x_polyarchy +
e_migdppcln + v2pepwrsoc +
(1 | year) + (1 + cohortmatch5_15 | iso3n) ,
data = data_imp1,
thin= n.thin,
chains= n.chains,
iter = n.iter,
warmup = n.burnin,
control = list(adapt_delta = 0.95),
cores = n.cores,
save_all_pars = TRUE)
print(m3.ols, digits=4)
ptm2 <- proc.time()
print((ptm2[3] - ptm[3])/60) # Stop the clock
install.packages("brms")
library(brms)
Sys.time()
ptm <- proc.time() # Start the clock
m3.ols <- brm(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
autocracy_5 + v2pepwrsoc_5 + autocracy_5*v2pepwrsoc_5 + v2x_polyarchy +
e_migdppcln + v2pepwrsoc +
(1 | year) + (1 + cohortmatch5_15 | iso3n) ,
data = data_imp1,
thin= n.thin,
chains= n.chains,
iter = n.iter,
warmup = n.burnin,
control = list(adapt_delta = 0.95),
cores = n.cores,
save_all_pars = TRUE)
print(m3.ols, digits=4)
save.image(file = "m3.ols.RData")
ptm2 <- proc.time()
print((ptm2[3] - ptm[3])/60) # Stop the clock
####################################################
# Macro parameters = 750 effective samples
####################################################
n.thin <- 2
n.chains <- 3
n.iter <- 1000
n.burnin <- 500
n.cores <- 3 # Use only this setting if you can spare multiple cores
#n.cores <- 1 # Local
print(n.thin)
print(n.chains)
print(n.iter)
print(n.burnin)
print(n.cores)
m3.ols <- brm(income_equality ~ age + as.factor(sex) + as.factor(education) +
as.factor(income_quintiles) + as.factor(unemployed) + as.factor(children) +
autocracy_5 + v2pepwrsoc_5 + autocracy_5*v2pepwrsoc_5 + v2x_polyarchy +
e_migdppcln + v2pepwrsoc +
(1 | year) + (1 + cohortmatch5_15 | iso3n) ,
data = data_imp1,
thin= n.thin,
chains= n.chains,
iter = n.iter,
warmup = n.burnin,
control = list(adapt_delta = 0.95),
cores = n.cores,
save_all_pars = TRUE)
print(m3.ols, digits=4)
ptm2 <- proc.time()
print((ptm2[3] - ptm[3])/60) # Stop the clock
#### "Gender Inequality and Authoritarian Regimes" ####
# authors: "Lars Pelke"
# date: 2019-09-16
# clear workspace
rm(list=ls())
#libraries
library(tidyverse)
library(ggpubr)
library(panelr)
library(readstata13)
library(texreg)
library(dotwhisker)
library(broom.mixed)
# set working directory
setwd("M:/Promotion/Paper/Gender Inequality/calculations/data")
#### Import datasets ####
##Import VDEM DATA
gender_data <- read_csv2("gender_data.csv")
##Resclaing Fuel Variables
gender_data$e_total_fuel_income_pcln <- log(gender_data$e_total_fuel_income_pc + 1) # plus 1, because zero values in original data
gender_data <- gender_data %>%
mutate(v2lgfemleg2 = ifelse(is.na(v2lgfemleg), 0, v2lgfemleg))
#### descriptive evidence ####
summary(gender_data$v2pepwrgen) # power distributed by gender
summary(gender_data$v2lgfemleg) # female legislators, set to missing when v2lgbicam is 0.
summary(gender_data$v2lgfemleg2) # transformed variable.
summary(gender_data$v2x_genpp) # women political participation (v2lgfemleg + v2pepwrgen)
gender_data <- gender_data %>%
mutate(v2lgqugen = ifelse(v2lgqugen>0, 1, 0))
summary(gender_data$v2lgqugen) # gender quota
gender_data$v2lgfemleg <- gender_data$v2lgfemleg/100
gender_data$v2lgfemleg2 <- gender_data$v2lgfemleg2/100
gender_data <- gender_data %>%
mutate(e_wb_pop_ln = log10(e_wb_pop))
#### Generate Sample ####
gender_data <- gender_data %>%
drop_na(v2lgfemleg2)
gender_data <- distinct(gender_data, country_id, year, .keep_all= TRUE)
gender_data2 <- gender_data %>%
drop_na(v2x_regime, year, v2pepwrsoc, e_wb_pop_ln, e_migdppcln, v2lgqugen)
gender_data3 <- gender_data2 %>%
drop_na(gwf_party, gwf_personal, gwf_monarchy, gwf_military)
#### Computing group_mean and de-meaned variables ####
library(sjmisc)
gender_data2 <- sjmisc::de_mean(gender_data2, v2x_regime, v2pepwrsoc, e_wb_pop_ln, e_migdppcln, v2lgqugen,
v2exl_legitideolcr_1, v2exl_legitideolcr_2, v2exl_legitideolcr_3, v2exl_legitideolcr_4,
v2exl_legitideolcr_0, v2exl_legitideol, v2exl_legitratio, v2exl_legitperf,
grp = country_id)
gender_data3 <- sjmisc::de_mean(gender_data3, v2x_regime, v2pepwrsoc, e_wb_pop_ln, e_migdppcln, v2lgqugen,
v2exl_legitideolcr_1, v2exl_legitideolcr_2, v2exl_legitideolcr_3, v2exl_legitideolcr_4,
v2exl_legitideolcr_0, v2exl_legitideol, v2exl_legitratio, v2exl_legitperf,
grp = country_id)
####################################################
# Bayesian Multilevel Models
####################################################
n.thin <- 2
n.chains <- 3
n.iter <- 500
n.burnin <- 200
n.cores <- 3 # Use only this setting if you can spare multiple cores
#n.cores <- 1 # Local
print(n.thin)
print(n.chains)
print(n.iter)
print(n.burnin)
print(n.cores)
library(brms)
## OLS
Sys.time()
ptm <- proc.time() # Start the clock
m1.ols <- brm(v2lgfemleg2 ~ year + v2exl_legitideolcr_1_dm +
v2exl_legitideolcr_1_gm + v2exl_legitperf_dm + v2exl_legitperf_gm + v2exl_legitratio_dm +
v2exl_legitratio_gm +
v2x_regime + (1 + year | country_id) +
(0 + v2exl_legitideolcr_1_dm + v2exl_legitperf_dm | country_id),
data = gender_data2,
thin= n.thin,
chains= n.chains,
iter = n.iter,
warmup = n.burnin,
control = list(adapt_delta = 0.95),
cores = n.cores,
save_all_pars = TRUE)
print(m1.ols, digits=4)
ptm2 <- proc.time()
print((ptm2[3] - ptm[3])/60) # Stop the clock
####################################################
# Bayesian Multilevel Models
####################################################
n.thin <- 2
n.chains <- 3
n.iter <- 1000
n.burnin <- 500
n.cores <- 3 # Use only this setting if you can spare multiple cores
#n.cores <- 1 # Local
print(n.thin)
print(n.chains)
print(n.iter)
print(n.burnin)
print(n.cores)
library(brms)
## OLS
Sys.time()
ptm <- proc.time() # Start the clock
m1.ols <- brm(v2lgfemleg2 ~ year + v2exl_legitideolcr_1_dm +
v2exl_legitideolcr_1_gm + v2exl_legitperf_dm + v2exl_legitperf_gm + v2exl_legitratio_dm +
v2exl_legitratio_gm +
v2x_regime + (1 + year | country_id) +
(0 + v2exl_legitideolcr_1_dm + v2exl_legitperf_dm | country_id),
data = gender_data2,
thin= n.thin,
chains= n.chains,
iter = n.iter,
warmup = n.burnin,
control = list(adapt_delta = 0.95),
cores = n.cores,
save_all_pars = TRUE)
print(m1.ols, digits=4)
ptm2 <- proc.time()
print((ptm2[3] - ptm[3])/60) # Stop the clock
